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ABSTRACT 

Orbital energy relaxation around a massive black hole (MBH) plays a key role in establishing the 
stellar dynamical state of galactic nuclei, and the nature of close stellar interactions with the MBH. 
The standard description of this process as diffusion in phase space provides a perturbative 2nd-order 
solution in the weak two-body interaction limit. We carry out a suite of A-body simulations, and 
find that this solution fails to describe the non-Gaussian short timescale evolution of the energy, 
which is strongly influenced by extreme events (a "heavy-tailed" distribution with diverging moments) 
even in the weak limit, and is thus difficult to characterize and measure reliably. We address this 
problem by deriving a non-perturbative solution for energy relaxation as an anomalous diffusion 
process driven by two-body interactions, and by developing a robust estimation technique to measure 
it in A-body simulations. These make it possible to analyze and fully model our numerical results, 
and thus empirically validate in detail, for the first time, this theoretical framework for describing 
energy relaxation around an MBH on all timescales. We derive the relation between the energy 
diffusion time, tg, and the time for a small density perturbation to return to steady state, t r , in a 
relaxed, single mass n(r) oc r~ 7 / 4 cusp around a MBH. We constrain the modest contribution from 
strong stellar encounters, and measure with high precision that of the weakest encounters, thereby 
determining the value of the Coulomb logarithm, and providing a robust analytical estimate for in 
a finite nuclear stellar cusp. We find that t r ~ 10ie ~ (5/32)Q 2 Ph/Nh log Q, where Q = M,/M+ is 
the MBH to star mass ratio, the orbital period Ph and number of stars Nh are evaluated at the energy 
scale corresponding to the MBH's sphere of influence, Eh — c^,, where is the velocity dispersion 
far from the MBH. We conclude, scaling Coo by the observed cosmic M./tJoo correlation, that stellar 
cusps around lower-mass MBHs (A/. < 10 7 M Q ), which have evolved passively over a Hubble time, 
should be dynamically relaxed. We briefly consider the implications of anomalous energy diffusion for 
orbital perturbations of stars observed very near the Galactic MBH. 

Subject headings: Galaxies: nuclei — Stars: kinematics and dynamics — Black hole physics 



1. INTRODUCTION 
1.1. Background 

Observations of our own Galactic Center (GC) re- 
veal that stars there move on Keplerian orbits in the 
potential of a central compact object, which is consis- 
tent with a massive black hole (MBH) of mass M. ~ 
4 x 10 6 M^ (lEisenhauer et al.l l2005t iGhez et al.1 [2TJ05I . 
l2008tlGillessen et al.H2009aU bD. It is widely assumed that 
the Milky Way is not unique in this sense, and that there 
is an MBH in the center of most, if not all, galaxies. 

The GC, with its relatively low-mass MBH, represents 
a sub-class of galactic nuclei where the dynamical relax- 
ation time is expected to be shorter than the Hubble 
time (Alexander 2007), so that the stellar distribution 
near the MBH could have had time to reach a relaxed 
state that is independent of initial conditions. More- 
over, because of its proximity, it is the only system to 
date where the dynamical state can be directly probed 
by observations. However, recent star counts in the GC 
(jBuchholz et al.|[2009t IDo et al.ll2009l ; iBartko et al1[20Tot) 
indicate that the radial density profile of the spectro- 
scopically identified low-mass red giants on the ~ 0.5 pc 
scale, thought to trace the long-lived population there, 
rises inward much less steeply (or even decreases) than 
is expected for a dynamically relaxed stellar population 
(IBahcall k Wolfl[il)76l [19771: 1 Alexander fc Hopmar] [2009l: 
iPreto fc Amaro-Seoane!l2010D . The simplest interpreta- 



tion is that, contrary to expectations, the GC is not dy- 
namically relaxed by two-body interactions. This could 
be due, for example, to a cosmologically recent major 
merger event with another MBH, whic h carved out a 
central cavity in the stellar distribution (Merritt 2010), 
or due to the presence of faster compe ting dynamical 
proce sses that drain stars into the MBH ( Madigan et al. 

MM)- 

This interpretation is however further complicated by 
a recent study of the diffuse light density distribution in 
the GC (|Yusef-Zadeh et al.l [2012h . which is thought to 
track the low-mass, long-lived main-sequence stars that 
are too faint and numerous to be individually resolved, 
and consist of the bulk of the stellar mass. In contrast 
with the stellar number counts, the diffuse light distri- 
bution appears to rise inward rather steeply down to the 
~ 0(0.01) pc scale, where destructive stellar collisions 
(not accounted for in the purely gravitational treatment 
of two-body relaxa tion), can plausi bly suppress the stel- 
lar density profile (| Alexanderl 1 1 9991 ) . 

These uncertainties about the dynamical state of the 
GC, and by implication, the state of galactic nuclei in 
general, have far-reaching ramifications. The possibility 
that the GC is unrelaxed severely limits the validity of 
extrapolating results derived from the unique observa- 
tions of the GC to other galactic nuclei. While a relaxed 
system can be understood and modeled from first prin- 
ciples, independently of initial conditions, an unrelaxed 
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one reflects its particular formation history, which is ex- 
pected to vary substantially from galaxy to galaxy. This, 
and the possibility that the GC lacks the predicted high- 
density central stellar cusp, have in particular crucial im- 
plications for attempts to understand the dynamics of 
extra-galactic gravitational wave (GW) sources and to 
predict their rates, since the low-mass MBH in the GC 
has long been considered the archetype of low-frequency 
extra-galac tic targets for planned spa ce-borne GW ob- 
servatories (jAmaro-Seoane et al.ll2007h . 

The study of dynamical relaxation around a 
MBH has a l ong hist o ry, b eg inning with the semi- 
nal stu dies of IPeeblesI (Il97l: IBahcall fc Wolfj (H976L 
19771 ): ICohn fc Kulsrudl (|1978D : IShapiro fc Marchantl 
who used various simplifying approximations 



Nelson fc Tremainel 11999ft to enable analytical and 
numerical approaches, and continuing more recently 
with studies of the time to approach steady state 
near a MBH in large A~-body simulatio ns (e.g. 
iPreto et al.ll2004HPreto fc Amaro-Seoandl2010D . A com- 
plementary approach of extracting the diffusion co- 
efficients (DCs) from stellar energy evolution in A- 
body simulations was carried out in several stud 
ies (E 
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remaind 119801: iRauch fc Tremaind 119961 : 



lEilon et all 120091 : iMadigan et alj 120111 ). leading to 
widely differing predictions about the dynamical state 
of galactic nuclei, es pecially very close to the MBH 
(jMadigan et al.l 1201 If ) . These open questions, and the 
fundamental importance of the process of orbital energy 
relaxation, motivate a re-examination of the analytic de- 
scription of energy relaxation, coupled to detailed com- 
parisons with carefully controlled simulations. 

1.2. Objectives and overview of this study 

In this study we use analytical considerations and high- 
accuracy Af-body simulations to address the following 
open questions: What is the relation between local en- 
ergy relaxation and the approach to a global steady 
state? What are the properties of evolution due to en- 
ergy relaxation on the short and intermediate timescales? 
How rare are extreme energy exchange events, and what 
is their contribution to relaxation? What is the physical 
meaning, and the correct value, of the Coulomb cutoffs 
introduced to control the formal divergences in the re- 
laxation rates? In practical terms, these questions all re- 
duce to the problem of relating the analytical / statistical 
description of local two-body interactions, to the global 
measurable properties of a stellar system. These issues 
must be addressed in order to better understand dynam- 
ical processes in galactic nuclei, improve their modeling, 
and reliably apply results from A^-body simulations to 
real physical systems. 

Our a pproach to these p roblems is complementary to 
that of IPreto et all (|2004f ). IPreto et all (2004T ) simu- 
lated a very large system with a non-Keplerian, and 
non-self-similar initial DF and potential and a rela- 
tively low mass ratio Q = M,/M+ over a relaxation 
timescale. IPreto et al.l (|2004l ) show that there is an agree- 
ment between the A^-body and Fokker-Planck (FP) sim- 
ulations by comparing the evolution of a small inner sub- 
system where the potential is dominated by that of the 
MBH and where the analytical/ numerical solution of 
IBahcall fc Wo!3 (1976ft (henceforth the BW cusp) is rel- 
evant. 



Here, we carry out a suite of small A-body simulations 
of a system with realistically high values of Q and a po- 
tential that is dominated by the MBH, designed specif- 
ically to simplify the analytic treatment and to allow 
complete modeling of the time-dependent evolution on 
short timescales. This allows us to rigorously test and 
verify the theoretical predictions, and in particular the 
scaling of the relaxation with energy and with the slope 
of the cusp. 

The choice of focusing on short timescales is to a 
large extent forced by the high computational cost of the 
simulations, since we model a high mass ratio system 
(Q = 10 6 , comparable to the mass ratio in the GC) with 
A~ > 10 4 particles which are strongly bound to the MBH. 
The short timescale behavior is however also of interest 
in itself, since dynamical processes near a MBH can be 
limited by the lifetime of the stars (for example, the mas- 
sive stars in the GC), or by fast destructive processes (for 
example, collisional destruction or tidal disruption). 

The elementary description of energy relaxation is by 
the energy transition probability, Ke,j (A£7) (Eq. I24p . 
which is the orbit-averaged probability per unit time, 
per unit energy, for the energy of a star to change from 
E to E + AE in a single scattering event. The resulting 
distribution of the accumulated change in energy over 
a finite time t, A t E, is represented by the propagator, 
WE,j(A t E, t), which describes the energy distribution of 
an initially mono-energetic system of test stars, after a 
time t. 

Past standard treatments of energy relaxation in stel- 
lar systems (e.g. IBahcall fc Wolfj|1976l) with the FP for- 
malism implicitly assumed that only the two lowest mo- 
ments of the energy transition probability (i.e. the 1st 
and 2nd DCs) play part in the evolution of the system. 
This is equivalent to assuming that the free propaga- 
tor is a Gaussian with width proportional to that 
is, energy evolves by normal diffusion. However, as we 
show, t he ~ 1/|A_E| 3 di vergence of the transition prob- 
ability (jGoodmanl I1983D implies that higher order DCs 
do play a role, which is dominant on short timescales, 
and diminishes with time, but still remains non-negligible 
even when the system relaxes to steady state. The impli- 
cations are that energy relaxation progresses somewhat 
faster than \fi, as an anomalous diffusion process. 

Our approach of measuring energy relaxation in A- 
body simulations on short timescales forces us to develop 
a robust statistical model for describing these timescales. 
By treating the problem in Fourier space, we are able 
to derive the full propagator explicitly. It can then be 
evaluated either numerically in non-perturbative form, 
or analytically in a perturbative manner. This enables 
us to show that the propagator is a heavy-tailed DF that 
evolves as ~ y/tlog(t) on timescales much shorter than 
the FP energy diffusion timescale. This short timescale 
behavior is an example o f an anomalous diffusion process 
(Metzler fc Klafteril200tift . One consequence is the much 
higher probability of events that would be expected to 
be extremely rare for a Gaussian probability distribu- 
tion. Another implication of the very slow convergence 
of this propagator to the central limit is that attempts 
to measure relaxation in relatively short A^-body exper- 
iments can yield very misleading results, which reflect 
more the method used, than the quantity probed. By 
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analyzing the scaling properties of the energy propagator 
we suggest a robust statistical estimator, using fractional 
moments, which is independent of the high energy cut- 
off. With it we measure the rate of energy relaxation on 
short timescales and confirm our analytical treatment, 
and in particular, the validity of the propagator and the 
DCs used in the standard FP approximation. Finally, we 
use the propagator to evolve a system by MC simulations 
up to steady state, and thereby confirm the approximate 
validity of the FP results on long timescales. 

The paper is organized as follows: In Section[2]we sum- 
marize the process of energy relaxation around a massive 
object in terms of a normal diffusion process described 
by the FP equation, and derive an analytical solution 
for the evolution of a stellar cusp to steady state. In 
Section [3] we present an elementary description of en- 
ergy relaxation in terms of the transition probability and 
its integrated propagator, and describe the emergence of 
anomalous diffusion that ultimately approaches normal 
diffusion on long timescales. In section @] we formulate 
a method to analyze our A^-body simulations and use it 
to confirm and calibrate our analytical description of re- 
laxation. In section [5] we briefly discuss the observable 
effects of anomalous diffusion on short period stars or- 
biting the Galactic MBH, perturbed by a dark cusp of 
stellar remnants. We summarize and discuss our results 
in section [6l 

2. ENERGY RELAXATION AS A DIFFUSION 
PROCESS 

The standard approach (jChand rasekhar 1943]) of treat- 
ing relaxation in a self-gravita ting system involves severa l 
simplifying assumptions fe.g. lNelson fe Tremaind ri999). 
The 2-body perturbations experienced by a test star are 
assumed to be independent of each other, that is, non- 
coherent. The encounter is assumed to be between two 
unbound stars, approaching each other from infinity. The 
perturbations are assumed to be local in two respects. (1) 
The exchange of energy and angular momentum occurs 
over a very short time at the point of closest approach. 
(2) The local conditions (stellar density and velocity) at 
the position of the test star, at the time of interaction, 
apply to all perturbers irrespective of their impact pa- 
rameter, which is equivalent to assuming an infinite ho- 
mogeneous system around the point of interaction. These 
assumptions cannot be satisfied for arbitrarily large im- 
pact parameters, and therefore require the introduction 
of a large-distance cutoff, which enters in the Coulomb 
factor A. Note however that it is far from obvious that 
these assumptions hold near a MBH. The potential there 
is dominated by the central mass, and on intermediate 
scales (close to the MBH, but far enough to ignore rela- 
tivistic effects) is approximately Keplerian, <f> — GM,/r; 
the stellar density decreases steeply with radius; and 
there is strong energy- dependent stratification of dynam- 
ical timescales. We discuss here the implications of these 
factors on energy relaxation near a MBH. Finally, it is 
commonly assumed that this Markovian process is also 
a diffusive processes, and can therefore described by the 
FP equation. 

2.1. The FP description of relaxation around an MBH 

Determining the stellar distribution in the vicinity of 
an MBH is a classical problem in stellar dynamics, first 



discussed in the context of an MBH embedded in a glob- 
ular cluster tiPeeblesl[l97i IBahcall fe Woljfl97(l . The 
general time evolution of the phase-space distribution 
function (DF) of test particles with mass m<r, /t(v, x, £), 
due to two-body interaction with field stars in a local 
velocity-mass dis tribution ^(v^ mp), is given by the 
master equation (|Goodmanlll985D 

-^/ T (v,x,i) = J dV/ T (v',x,t)X(v',v) 



■/ T (v,x,i)/dV%v'), (1) 



where D/Dt = (d/dt)+v-(d/dx)-(dcf>/dx)-(d/dv) and 
AT(v,v') is the transition probability per unit time per 
unit velocity volume that a test star changed its velocity 
from v to v' = v + Av as a result of a single encounter 
with a field star 



4G 2 f f 

i4: ( v ' v ') = TA — is / dm Fm 2 F / 

|Au| a J J vr . Av 



f F w + v + A- 



UlT + TTlF 

T 

2m f 



,m F 



(2) 



This description of stellar dynamics as a relaxation pro- 
cess assumes that changes in energy and angular momen- 
tum result from uncorrelated velocit y deflections, which 
are local in time and space (e.g. iNelson fc Tremaind 
1999; iFreitad [2008T ). Dynamical mechanisms that in- 
volve long timescale correl ations, e.g. resonant relaxation 
(|Rauch fc Tremaind [l996h . are not taken into account 
here. 

Assume that fr and fp are functions only of orbital 
energy E and time t\ the potential is dominated by a 
massive central object <\> = GM m /r; and the field stars 
and the test stars have the same mass, M t . The master 
equation ([1]) is then reduced to an equation of energy 
only (here and below the energy and potential of bound 
stars is defined as positive), 

dN( f t ,t] = J dAE{N(E - AE)K E _ AE (AE) 

-N(E)K E (AE)}, (3) 

where K E ( Ae) is the orbit- and eccentricity-averaged en- 
ergy transition probability, and N(E) is the stellar num- 
ber density in energy space 



N(E,t)=4Tr 2 J^E)P(E)f(E,t). 



(4) 



By assuming weak encounters, AE <C E, and expand- 
ing N (E — AE) and K E _ AE {AE) around E — AE = E, 
the master equation can be express ed as a Kramers- 
Moyal expansion (e.g. Weinberg 20011) 

dN(E,t) ^ 1 / d \ n A 

where the DCs are the moments of the transition proba- 
bility K E (AE), 

,AB max 

r£(£)= / AE n K E (AE)dAE, (6) 

JA_E min 

and where A = AE maK / AE min . For n — 1,2, P^ is 
proportional to the Coulomb logarithm log A (Eqs. [Hjl 
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1 1 1 p ; for higher moments the dependence is much weaker, 
oc (1 — A - L("— !)/ 2 J ) (Eq. [3i]) j and can usually be ignored. 

By truncating the sum after the first two terms, the 
Kra mers- Moyal expansion is reduced to the FP equation 
(e.g. iRosenbiuth et al.lll957ft 



dN(E,t) 1 d 2 



dt 



2dE 2 
d 



{N(E,t)T^E)} 



dE 



{N(E, t)Tf(E)} = 



8F(E) 
dE 



(7) 



where F(E) is the stellar flux. 

The FP approximation is valid only if the higher or- 
der terms can be neglected. This is usually justified 
by arguing that the n > 2 terms of the expansion are 
of the order of the n > 2 DCs themselves, which are 
smaller than the first two DCs by a fact or of 1/ log A (e.g. 
Binney & Tremaine 2008; Henon 1960). However, as we 
show in SectionGOJ such a comparison is meaningful only 
when d/dE ~ 1/AJ5 ~ l/E (Eq. |5TJ), i.e. the change in 
energy has accumulated to an order unity change. This 
is a valid assumption only over long timescales; on short 
timescales the system cannot be described by the first 
two DCs only. 

It is further commonly assumed that stars that are 
unbound to the MBH (E < 0) are drawn from a Maxwell- 
Boltzmann distribution, 



/(£) = 



^3/2 



E < 0, 



(8) 



Mr 

where no and <7o are the stellar number density and veloc- 
ity dispersion of this simplified model of the host galaxy. 
The inner energy boundary is determined by the max- 
imal orbital energy E& where stars can survive against 
any of the various disruption mechanisms that operate in 
the extreme environment near a MBH (e.g. tidal disrup- 
tion, tidal heating, orbital decay by gravitational wave 
emission, stellar collisions), 



f{E) = 0, E < E d 



(9) 



The steady state solution {dN/dt = 0) of the FP equa- 
tion can be simply o btained if a cu sp solution f(E) — 
CE p is assumed (e.g. iPeeblesI [l97l . Although this so- 
lution does not satisfy the boundary conditions (Eqs. [5J 
[5]), it is in fact a reasonable approximation of the full 
solution far from the boundaries at -C E <C Ed, (e.g. 
Bahcall fe Wolf I1976T) . For a full analytic solution see 



Keshet et al.l (|2009ft 



Approximating the DF by a pure power-law, — 1 < p < 
1/2, the DCs are 



4(1 -Ap) \ogAN{< E) 



l+p 



Q 2 P{E) 



E oc E pJl 



and 



16 (3 -2p) log A A^(< E) 



(l+p)(l-2p) Q 2 P(E) 
The steady state solution is then (Eq. [7]) 

2 



E 2 oc E pA 



(4p-l)(4p-3)- 



2p 



0, 



(10) 



(11) 



(12) 



whose solutions are either p — 3/4 (Peebles 1971), which 
is unphysical since T^(i?) diverges with Ed — > oo leading 



to hu ge stellar flux away from the MBH (|Bahcall fe W olf 
or p = 1/4, the so-called BW cusp solution 
(|Bahcall fe Wolfl 11976ft . The BW cusp solution has 
zero net flow (F = 0) since the term in the equa- 
tion that contains the second DC is independent of en- 
ergy, and since there is no drift (jLightman fe Shapiro! 
11977ft . This solution was confirme d by integrating the FP 
equation, both numerically fe.g. iBahcall fc Wolj 119761 ; 
iCohn fc Kulsrudl 119781: iPreto et al.l 12004ft and statisti- 
cally (e.g. iShapiro fc Marchantl 11978ft . bv Monte-Carlo 
(MC) simulations. 

2.2. Relaxation to steady state 

Stellar systems out of equilibrium evolve by redis- 
tributing their energy and angular momentum. When 
this proceeds by diffusion, the timescales associated with 
the various DCs generally reflect the dynamical evolution 
timescale. One such timescale is the energy diffusion time 

t E {E) = E 2 /T${E) (13) 

which expresses the timescale for a test star to change 
its orbital energy by order unity. Another commonly 
used quantity is the Chandrasek har relaxation time 
(|Chandrasekhanll942USpitzerlll987ft . 



;(a«h)s 



0.34 



=V3cr 



G 2 r 



■An A 



(14) 



where the second equality expresses the timescale for 
a typical test star with velocity v = V3a, to change 
its kinetic energy by order unity, in a Maxwellian, 
isotropic and homogeneous stellar system with stellar 
mass m, number density n and velocity dispersion a . The 
Coulomb logarithm is defined by the ratio of maximal 
and minimal impact parameters In A = ln(p max /p m i n ), 
where p mlu = Gm/a 2 (small angle approximation), and 
Pmax is the size of the system. 

These two timescales are related, but different, since 
tch is locaQ in physical space (r, v), whereas is local 
in (E, J) space, and hence is orbit-averaged over all (r, v) 
values along the orbit, and often samples a wide range 
of conditions in the stellar cluster. It is therefore to be 
expected that the values of the two timescales can be 
quite different, even though they express the same basic 
property of the system: systems with short relaxation 
timescales evolve rapidly, and those with long timescales 
evolve slowly. For example, when the two timescale are 
compared at E ~ er 2 , typically ich ^> (Eq. [21]) . 

The diffusive timescales cannot be used in themselves 
to estimate how long it takes an out-of-equilibrium sys- 
tem to approach steady state, assuming that such ex- 
ists. This generally depends on the specific initial con- 
ditions, and the operative definition of convergence to 
steady state. Numeric experiments show that in many 
cases tch, evaluated at the radius of MBH influence, is 
a reasonable estimator for the relaxation time to steady 
state, ty . 

We now proceed to use the FP equation to demonstrate 
analytically why that is so, and to derive the relation be- 
tween ie, tch and t r in a power- law stellar cusp around 

1 The distinction between local and global is irrelevant for tc^ 
when an infinite homogeneous system is assumed. 
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a MBH. In Section |3~31 we verify that our analysis is not 
limited by the approximate nature of the FP treatment, 
by numerically integrating the system using MC simu- 
lations where the individual energy steps are generated 
using the isotropic propagator We (At £7). Moreover, we 
show that strong encounters, which are ignored in the 
FP treatment, reduce the time to reach steady state by 
a factor of only < 4/3. 

Consider a background stellar population in 
steady state with a power- law cusp Nb g (E) = 
(5/4)N(E/E p )- g / 4 , where N is the total number 
of stars with E > E h . The DCs (Eqs. [lOj HTJ are 
r£(E) = and T%(E) oc E 9 / 4 . Assume a small initial 
perturbation superimposed at energy Eq on the steady 
state background, N P (E) = N±5(E — Eq). This system 
describes, for example, test stars that are scattered by 
a steady state distribution of field stars. Over time, 
the perturbation will spread out until it relaxes to the 
oc _E -9 / 4 steady state distribution. 

The dimensionless FP equation for the evolution of N p 

is 

where x = E/Eh and r = t/tE{Eh)- We assume that 
the number of stars with E > Eh is constant (i.e. there 
is a reflective bou ndary at Eh). The perturbation then 
evolves as (see Eq. ICllj) 

oo 

V (x,r)= ^$-1= X)Q«(aOQn(*o)e-£» T/128 GL6) 

n=l 

where N™{x) = N bg {x)/N, Q n {x) 

V'^ J 4(j5,n^~ 1/8 )M(.?5 : n), J n {x) is the Bessel 
function of the first kind and j n ,m is the m'th zero of 
J n (x). Similar solutions for a power-law background 
cusps with p > are presented in Appendix [C] 

The time it takes for the small perturbation to decay, 
and for the system to converge back to its steady state, 
depends on the location of the initial perturbation, and 
on the operative definition of the convergence. The dy- 
namical evolution back to steady state is slower for higher 
Xq (closer to the MBH), in spite of the fact that the re- 
laxation time decreases in that limit as 4e(:e) ~ x~ p 
(p — 1/4). This counter-intuitive result can be un- 
derstood by considering a perturbation of mass AM at 
x > xq, where the enclosed steady state mass scales as 
Xq 3 ^ 2 . It then follows the fractional perturbation scales 
as 8M oc x^ 2 p , and the rate for relaxing the pertur- 
bation scales as SM/tE oc x^ 2 . We therefore conserva- 
tively define the relaxation time as the time it takes a 
5-function perturbation in the limit xq — » oo to decay to 
an order unity perturbation, i](xo.t r ) = 1, where 

oo 

„(*)= Urn £Q2(zo)e-^V/(i28 tB ) 

n—1 

1 J tl »-jg,it/(128te) (i 7 \ 

1128r(5)2 J|(i 5 ,i) ' [ ' 

where here T is the Gamma function. This corresponds 



to the relaxation time 

tr = |l log (ii28r(5) 2 j/So) tE * lhltE - (18) 

The conclusion that t r <~ O(10)£e is robust, since it 
mostly reflects the time it takes the perturbation to prop- 
agate down to x = 1, rather than the exact details of the 
convergence criterion. For example, an alternative defini- 
tion of t r in terms of evaporation, leads to a similar result. 
Consider the evaporation of a perturbation through the 
outer boundary x = 1, and define the relaxation time to 
steady state as the time it takes the number of stars still 
in the system to decrease by a factor P+ = 0.1 (the con- 
clusions depend only logarithmically on the exact value 
of P+). In this setting the probabil ity to find the test star 
in the system is given by (see Eq. lCl5[) 

N p (x > 1,T) _ _2_ J 5 (j4,n)- / 4(j4,"/ a: /8 ) c -7, 2 „r/128 
N * V*0 ~[ k.nJlijA.n) 

( 19 ) 

and the relaxation time (in the limit x — > oo) is given 
by 

3i,i V8J 5 (j 4 ,i)r(5) J 

similar to the relaxation time we obtained before. 

This simple model can also explain the behavior 
seen in a variety of dynamical simu l ations of relax- 
ation around a MBH (IBahcall fc Woli 111 976ft Figure 2; 
iHopman fc Alexanderi (|2006f ) Fi gure 1: IMadigan et aT] 
(2011) Figure 30): the initial relaxation to steady 
state is extremely rapid (super-exponential), while 
at later times it slows down to exponential decay. 
Eq. (TIT)) indeed shows that the perturbation con- 
verges to steady state as a series of exponential de- 
caying terms with progressively faster rates 128/jf n = 
{0.601, 1.19, 1.93, 2.81, 3.86, . . . }. Thus, the higher 
terms vanish rapidly in a combined super-exponential 
rate, leaving last only the first and slowest term, which 
decays exponentially on a timescale > tE- 

The relaxation time to steady state, t r , can be related 
to the Chandrasekhar time ich in the host galaxy (ide- 
alized here as an isotropic and homogeneous system). 
Beyond the MBH radius of influence at rh — GM % ja\, 
where cto is the asymptotic velocity dispersion and the 
stars are unbound to the MBH (E < 0). The veloc- 
ity distribution is Maxwellian, the DF is given by Eq. 
J8|), and therefore tch is given by Eq. (fl%| . with the 
asymptotic stellar density no. The typical orbital energy 
at rh is ~ Cq/2. Closer to the MBH, within the ra- 
dius of influence where the dynamics are Keplerian, the 
steady state DF for bound stars with E > [~ rfc/2) 
is / w 2/(0)(£/cto) P (e.g. IBahcall fc WolJfl976h . The 
energy diffusion timescale is then (Eqs. [131 ITTj) 

It then follows that t r ~ 10t E ~ T ch /4 (Eqs. [U |2T|) . 
The large numeric pre-factors that relate these different 
diffusion timescales highlight the point discussed above. 
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While these timescales all generally reflect the tendency 
of the system to evolve, they measure different aspects 
of it: tch is a local quantity in physical space, estimated 
outside the evolving region; tE is an orbit-averaged quan- 
tity, estimated at some typical radius inside the evolving 
region; t r describes the timescale to reach global steady 
state. In many cases of interest, t r is not very different 
from the age of the stellar system (or the Hubble time, 
a commonly assumed upper limit). Inferences about 
the dynamical state of such systems should be mind- 
ful of the large differences between the meanings and 
values of the various diffusion timescales. The result 
tr ~ Tch/4 is consistent with numerical studies that show 
that a stellar system around a MB H approaches steady 
state in about T Ch /3 to T Ch / 2 (e.g. lBahcall fc Woljfl97a 
Hop man fc AlexanderlfeOOfl ) . 

The scaling of t r with MBH mass M, can be obtained 
from t he relation t r — lOtp;, the empirical Mia re- 
lation (iFerrarese fc Merrittl 120001 ; iGebhardt et "all 120001; 
ITremaine et all 120021 ), " 

M. = M (a/a ) a , (22) 

and the stellar number parametrization N (> E) = 

Q(E/(,a 2 y 5/4 , where £ is an order unity fac- 
tor. The time to steady state is then t r = 

(5V2tt/64) r 5/4 (M /M.) 3/a (GM./af) Q/logQ. For 
a = 4, M n = 1.3 x 10 8 M o , and a = 200 kms" 1 
(ITremaine et al.ll2002l ) 



io 5 - 



L. = 7.1 



1/4 



(O/io 7 ) 



5/4 



(logQ/loglO 7 ) 



(23) 

This indicates that stellar cusps around lower-mass 
MBHs (M. < 1O 7 M with M* = 1M Q ), which have 
evolved passively over a Hubble time, should be dynam- 
ically relaxed. 

3. BEYOND THE DIFFUSION APPROXIMATION 
3.1. Kinematic limits on energy exchange 

The transition probability encodes the full descrip- 
tion of energy evolution in the system due to 2-body 
encounters on all energy scales. For hyperbolic Keple- 
rian interactions between two stars, its general form is 
K r ,E ~ F r ^E(AE)/\AE\ 3 , where F r _E is a non-diverging 
function of AE, the energy change in an encounter (see 
Appendix |A"|) . The singularity of K r ^E at AE — > is 
never reached in practice, because the finite size of the 
system sets a lower limit on AE. In fact, an even stronger 
limit is implied by the restriction to hyperbolic orbits. 
A test star on an orbit with energy E has radial or- 
bital period P{E) ~ GM.E~ 3 I 2 and a typical veloc- 
ity v ~ \~E. An encounter with another star at im- 
pact parameter b leads to an exchange of energy of order 
AE ~ GAE/b and lasts a time ~ b/v. Since the weakest 
encounter that can still be approximated by a hyperbolic 
one, must have b/v < P, it follows that AE m j n ~ E/Q 
(i.e. Ae^jn - 1/Q for Ae = AE/E) (|Bahcall fc Wolfl 
1976). This approximate effective \ow-AE cutoff as- 
sumes that the contribution of distant and slow inter- 
actions is negligible. It is estimated below empirically by 
iV-body simulations (see Section |4~2|) . 




0.001 



0.0 0.5 1.0 

Fig. 1. — The transition probability K T: e (Eq. IA2[| for the case 
of a/r = 1 and a BW cusp (p = 1/4), as function of Ae. The full 
transition prob ability (solid line) is compared to the approximate 
form (Eq. IA3|l (dashed line). 

Kinematics in troduce an inherent high-A-E limit 
([Goodman! Il985h . Strong encounters occur in the im- 
pulsive limit, where the distance between the test star 
at r and the background star at r^ is very small. The 
energy is then exchanged over a short time and the test 
star's position remains almost fixed. In that limit, the 
energy one star can transfer to the other is at most 
its total kinetic energy. Let the test star be on an 
initial Keplerian orbit with semi-major axis a, energy 
E = GM,/2a eccentricity e and the eccentric anomaly 
oj, is at radius r = a (1 — e cos w) at the moment of the 
encounter. The most bound orbit it can be deflected 
to at fixed r has a' min — r/2 and E' max > E. Thus, 
the maximal change in the relative energy of the test 



star is Ae 



(+) _ 



mzx\E' - E\/E = (2o/r-l). Con- 
versely, the least bound orbit the test star can be de- 
flected to at fixed r involves kicking a background star 
with an initial energy Ef, to a maximally bound or- 
bit with a' b — rb/2 s» r/2. Since typically the least 
bound (or unbound) background stars have Eb ~ 0, 
the maximal change in the relative energy in this case 

is Ae,n a x w 2a/ r. Figure {1} shows the full transition 
probability for a specific value of a/r and a power-law 
DF; the kinematic high-AE 1 limit is clearly seen. The 
high-AE limit is ultimately set by physical (inelastic) 
collisions between stars of finite size, where orbital en- 
ergy is not conserved due to tidal interactions, mass loss, 
or stellar destruction. Orbit-averaging over the phase of 



the test star yields ( Ae 



(+) 

max. 



1 and ( Ae 



2, close 



to the commonly assumed value Ae max = 1. It is inter- 
esting to note that these averaged kinematic limits are 
consistent with the small angle regime, conventionally 
defined by an impact pa rameter b > 2GAE/ (v) 2 (e.g. 
iBinnev fc Tremaindl200ct ) . which corresponds to Ae < 1 

Note however, that in the limit e — > 1, Ae diverges, and 
therefore the eccentricity-averaged transition probability 
has finite contributions from arbitrarily large Ae. In par- 
ticular, the isotropic average transition probability, and 
the resulting isotropic propagator, involve integrations 
over Ae — » oo. 

Note that the actual value of Ae max may be lower 
than the kinematical one, when destructive physical col- 
lisions between stars play a role. In that case the max- 
imal energy exchange is given by Ae^ ax ~ Aa/QR* m 
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45(Q/4xl0 6 ) (a/pc) (-R*/-R ) \ where i?* is the 
stellar radius. However, this is relevant only extremely 
close to the MBH. For example, in the Galactic center 
Ae max < 1 for a < 22 mpc, and even then the effect of the 
collisional limit on the relaxation of the surviving stars 
is small; for example at a ~ 0.22 mpc, log A is reduced 
by only a factor of two. 

3.2. Relaxation on short timescales 

In order to measure short-timescale evolution in N- 
body simulation, it is necessary to derive the en- 
ergy propagator in the short time limit. Consider a 
test star with an initial energy E and angular mo- 
mentum J, where the energy and angular momentum 
change as a result of two-body interactions over a pe- 
riod of time P(E) < t < t E . In the small deflec- 
tion angle limit (equivalently, AE / E < 1), the or- 
bit averaged energy transition probability, K Ei j(AE) = 
[ dA JKp!_j( AE, A J), can be a pproximated as (See Ap- 
pendix |A] and [Goodman] ()1983h ) 



K E j(AE) = 



T 2 {E, J) + J)AE + Q(AE 2 
2|A£| 3 



(24) 

for AE min < \AE\ < E, where T 2 (E, J) and T 1 (E, J) are 
given by Eqs. (| A6|) and (|A7[) . Note that Ti^ are shape 
parameters of the transition probability, Ke,j (AE), and 
are not by themselves the DCs T^ 2 , which are related to 
them by 

r^ 2 (B,j) = r 1)2 (£; ) j)iogA. (25) 

This distinction is important because log A ~ log Q 
is typically a large, O(10) factor, which is determined 
by our assumptions that strong encounters are ignored. 
Later is section l3~3l we will address this assumption and 
estimate the contribution of strong encounters. However 
we show below that the short timescales evolution of the 
system does not depend on the high AE 1 - limit. 

Assume that t is short enough so that the back- 
ground remains fixed, and changes in E and J are small 
enough so that Ke,j(AE) is approximately constant 
along the trajectory of the test star in (E, J) space. Let 
WE,j(A t E, t) be the energy propagator. In that case the 
master equation is 



dW EJ (A t E,t) 
dt 



J K E j(AE)W E j(A t E~ AE.t)dAE 
-q 2b (E,J)W E! j(A t E,t), (26) 



where W Ei j(A t E, 0) = 5(A t E) and q 2b is the total rate 
at which 2-body encounters scatter the test star 

q 2b (E, J) — f dAEK E .j(AE) - l~£ 2 \ 2 (E, J) . 

(27) 

The energy propagator W E j(A t E, t) can then be used 
to integrate the master eq uation in a MC scheme (e.g. 
IShapiro fc Marchantlll978l ). 

The Fourier transform of Eq. (|2"6"|) is 

(^-K Ei j(k) + q 2b (E,J) \ W EJ (k 7 t) = 0, (28) 



where the solution is simply 

W E j(A t E,t)= ^ e ikA t E e (K E ,j(k)- g2b )t (2Q) 

J-oo 27r 



and K E j(k) is the Fourier transform of K E j(AE) 



K J ,E{k)=q 2b -ikT^- 1 ^-T$ 



Hfc) ? 

n— 3 



,(30) 



where are the n'th (non-central) moments of 

Kj t E(AE). The higher odd and even moments are 



1 2n+l 



E 2n T l 1 - A" 



2n 



Ae 2n 

max 



2n+2 



2 II 



Ae 2n 



(31) 

Note that Eq. (|30f has an analytical e xpre ssion in 
terms of cosine and sine integrals (see Eq. IB2|) . Thus 
a full, non-perturbative, evaluation of the propaga- 
tor We j(AE,t) is obtained by numerically evaluating 
Eq. (|29jl . This is equivalent to including in the master 
equation DCs of all orders (i.e. all moments of the transi- 
tion probability) , albeit calculated using an approximate 
transition probability, which incorporates the large angle 
limit only as an effective cutoff. This is to be contrasted 
with the standard FP treatment where only the first two 
DCs are included (often also only as approximations). 

Neglecting the contribution of the n > 2 moments in 
Eq. (|30|) results in the FP propagator 



WE P j(A t E,t) = 



1 



(A t E-r' v t)^ 



(32) 



which can be related to the non-pertubative energy prop- 
agator (Eq. HU) via Eq. J30]), 



W E ,j(A t E,t) = W^j(A t E,t){l 

1 1 1 ETi 
~23! logAT - 




where He n (x) are the Hermite polynomials and 
t E {E,J) = E 2 /(T 2 (E, J) log A) is the energy diffusion 
time. Therefore, for t > t^/logA, the higher moments 
can be neglected, as expected from the central limit the- 
orem. 

On short timescales, the deviation from the FP prop- 
agator can be significant, as seen in Figure [5] Figure [3] 
shows the evolution of the energy propagator. On short 
timescales (t < 0.01^), the FP approximation overesti- 
mates the contribution of small energy changes and un- 
derestimates the contribution of large energy changes to 
the accumulated change. Thus, formally rare events (i.e. 
> 3<T for a Gaussian) have a much larger probability than 
expected by the FP treatment. For example, the proba- 
bility of a 5<7 Gaussian event (2.9 x 10~ 7 ) can be actually 
as high as 1.6 x 10~ 2 . 

It should be emphasized that the non-Gaussian, super- 
diffusive short-timescale behavior occurs in spite of the 



Fig. 2. — Numerical evaluation of the propagator (Eg, I29H with AE n 



E and AE„ 



E/10 G for several time lags (solid lines). 



On the left the pdf is normalized by a 



it/tj5 as expected from the FP approximation (Eq. 1321 ) . on the right the pdf 



is normalized by a = ^JtT^it) as expected from the asymptotic analysis (first term in Eg. 1341 1. A normalized Gaussian A^(0, 1) (dashed 

line) is shown for comparison. The non-perturbative propagator converges to the approximate FP one by t > O.Ue- The convergence is 
faster in the central part of the distribution. The short timescale evolution of the central part is well approximated by a Gaussian with 

a = ytT^ (t), a clear demonstration of anomalous-diffusive evolution, which is faster than the a <x \/i of a standard random walk process. 




Fig. 3. — The probability P(> x) for a |A t _E| > xa FP event 
normalized by the probability to have such an event in a normal 

distribution, where a Fi 

tor. The FP propagator underestimates large a events, but over 
estimates the small a events. 



EJtT% denned by the FP propaga- 



fact that the physical energy transition probability (in 
contrast to its simplified power-law approximation) is in 
fact bounded by Ai? max , and therefore does have non- 
diverging moments. This seeming contradiction is rec- 
onciled by noting that on short timescales, this cutoff is 
irrelevant, since it is too far away in the wings of the 
distribution — the system does not evolve enough to "be 
aware" of its existence. This is a typical behavior for 
distributions whose variance would diverge were it not 
for p hysical cutoffs, for examp le truncated Levy-flights 
(e.g. iZumofen fc Klafterl 119941 ). Only on long enough 
timescales, when repeated scattering populates the wings 
all the way to the cutoffs, does the variance converge, the 
conditions for the application of the central limit theorem 
are satisfied, and the propagator converges to its asymp- 
totic a oc \fi behavior, independently of the functional 



form of K E j(AE). 

In Appendix |B] we derive the effective untruncated 
(Ae max — > oo) propagator W% j(A t E, t). Since it in- 
cludes arbitrarily strong energy exchanges, we refer to it 
at the "strong propagator" (SP). This propagator is given 
by a heavy tailed probability distribution 



W^j(A t E,t) = Wlj(A t E,t) 
whose leading term is a Gaussian 



Wlj(A t E,t) = 



Wlj(A t E,t) 
4 log L{q 2b t) 



V27rfr£(i) 



where 



rfaW = ri, 2 (£, J) logL(q 2b (E, J)t) 



(34) 



(35) 



(36) 



and L (x) w Ve 3 ~ 2 ~ fE x, where 7^ is the Euler constant. 
On long timescales approaching tg, L(q2bt) ~ A. The 
second, correction term in (|34[) is the higher order (in 
1/ logL (q2bt)) even function that contributes an asymp- 



totic tail of - \A t E\~ i as \A t E\ 



00, 



W} 



iA t E,t): 



A t g-trf(t) 



(37) 



where (x) is defined in Eq. (|B7|) . The next order odd 
correction function corresponds to the drift term (see Eq. 
IB12[) which is small on short timescales and is omitted 
here. Unlike a Gaussian, which scales self-similarly as 
o~{t) , the propagator in Eq. (|34[) is not self-similar and 
thus does not have a single scaling. However, its cen- 
tral part scales like ^/tT^it), which increases with time 
faster than \/i 



2 , that is, in a super-diffusive manner. 
On longer timescales, when the cutoffs become important 
and the limit A — ► 00 is not valid, Eq. (|3"4"|) is no longer 
a good approximation of the propagator. Note that on 
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all timescales this approximation is valid only for some 
central part of the propagator and it is never a good ap- 
proximation of the very far tails. Therefore this approxi- 
mated propagator can not be used in the MC simulation 
even if the time steps of the MC are very short. 

3.3. Relaxation on long timescales 

The steady state configuration of a system, if such is 
attained on cosmologically relevant timescales, has a spe- 
cial significance since it is independent of the initial dy- 
namical conditions, which are usually unknown, and is 
typically the most likely configuration to be observed. It 
is therefore important to understand how, and on what 
timescale, a system reaches steady state, and how short- 
timescale evolution by anomalous diffusion asymptoti- 
cally approaches normal diffusion, as described by the 
FP approximation. 

On short timescales, rare strong encounters (Ae > 1) 
can be neglected. However, as time progresses, even 
these rare events contribute significantly, and the kine- 
matic limit (Section 13. ip must be taken into account. A 
key issue is to estimate their effect on the relaxation rate. 
In principle, the exact propagator (EP) provides a full 
description of all encounters. However, from a practical 
point of view it is difficult to implement it in an analytical 
or numerica l scheme, due to the required triple integra- 
tion (but see Goodman 1983). Instead, as we show below, 
the EP can be bracketed by two limiting approximations. 
The first, the SP W^(A t E) (Eq. [34]), over-estimates the 
contribution of strong encounters by taking the weak- 
limit transition probability (Eq. I24p , and integrating it 
to Ae max — > oo, i.e. beyond the weak limit (Eq. IB2[) . 
The second, the weak propagator (WP) We{\E), in- 
tegrates only up to Ae max = 1, where the weak limit 
is valid (Eq. l30l) . As argued below (see Figures 01 [5| , 
the SP, WP and FP propagator all lead to evolution to 
the same steady state, on similar timescales. This then 
proves that so must the EP, and that the contribution of 
strong encounters to reaching steady state is not crucial. 
Table [1] summarizes the properties of the four isotropic 
propagators that enter our analysis of relaxation. 

The inclusion of strong encounters in the transition 
probability can only accelerate relaxation. Therefore, 
the SP is expected to lead to the fastest relaxation, fol- 
lowed by the EP, while the WP and the FP are slowest, 
since they assume the weak limit. The question of the 
existence and values of the moments of the propagators 
reflect the form assumed for the transition probability, 
and in turn are related to the nature of the long-term 
evolution of the system. The weak limit assumption in- 
troduces a cutoff on the energy exchange, which ensures 
that all moments exist. In particular, the neglect of all 
n > 2 moments assumed for the FP propagator, implies 
that it must be a Gaussian. Both the WP and the ex- 
act one are guaranteed by the Central Limit Theorem 
to converge to a Gaussian. The WP, by construction 
(Ae max = 1) converges to the FP Gaussian. This is not 
the case for the EP, which takes into account the actual 
kinematic limits. These can be shown to correspond in 
the asymptotic (Gaussian) long timescale limit to an ef- 
fective Ae^ x ^> 1. However, as we demonstrate below 
(Eq. 155]) . this limit is approached only on timescales 
much longer than the time for the stellar system to re- 
lax to steady state, and therefore has little effect on the 
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Fig. 4. — Time evolution of the perturbation's DF f p (x,t) for 
initial perturbation at xo = 100. The system is integrated by 
MC, with individual energy steps drawn from the FP propagator 
(triangles), the isotropic WP WE(AtE) (crosses), and the isotropic 
SP W§,(A t E) (asterisks). The analytically derived DF f p [x, t) (Eq. 
1161 is shown for comparison (solid lines) . All propagators converge 
to the FP DF for t > tg. The asymptotic steady state f p (x) ac X 1 / 4 
is also shown for reference (dash-dotted line). 

actual time to reach steady state. This even holds for 
the SP (Ae max — > oo), where all the moments diverge 
and the propagator never converges to a Gaussian with 
width oc \ft. 

The long-term evolution of the system can be de- 
scribed by repeated application of the propagator on 
short-enough timescales. Computationally, this can be 
realized by a MC simulation, which follows the evolution 
of the orbital energy of individual stars as it changes by 
small random increments, drawn from the distribution 
function described by the propagator. Figure 0] shows 
snapshots from a set of such simulations. In each, a dif- 
ferent propagator is used to follow the evolution of an ini- 
tial S- function distribution of test stars around a MBH, 
assuming a relaxed stellar background, and a reflecting 
boundary condition at the outer edge of the system. In 
all cases, the perturbation relaxes to the steady-state BW 
solution after t r ~ 10^. This provides numeric valida- 
tion of the analytical treatment of this process (Section 
12. 2|l . As expected, the evolving DFs resulting from the 
WP and SP are significantly broader at early times than 
the FP one. This reflects their heavy-tailed probability 
distributions (e.g. Figure [5]). 

Another demonstration of the near equivalence of the 
different propagators is shown in Figure [S] where the 
different propagators are used to follow the evolution of 
an initial 5-function distribution, assuming an absorb- 
ing boundary condition at the outer edge of the system. 
In this configuration, the time to relaxation can be de- 
fined as the time it takes the system to evaporate some 
large fraction of its stars. Since we are interested here 
mainly in comparing the evolution due to the different 
propagators, our conclusions do not depend on the ex- 
act fraction. For definitiveness, we choose 0.9. Both the 
FP and the WP drive evaporation on almost the same 
timescale, while the SP drives it only slightly faster, by 
a factor of - 1.1 (for Ae min = Q- 1 = 10" 6 ). 

This modest acceleration in the evaporation rate by the 
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TABLE 1 

Summary of properties of isotropic propagators 





Strong (SP) 


Exact (EP) 


Weak (WP) 


Fokker-Planck (FP) 


Strong encounters 


Overestimated 


Exact 


Not included 


Relaxation time 


MSP) < 


ir(EP) < 


i r (WP) = 


tr(FP) 


Moments 


All diverge 


Only n = 1, 2 finite 


All finite 


n = 1, 2 finite, n > 2 zero 


Gaussianity 


Only central part, a 2 ~ tlogt 


On long timescales, <r 2 <x t 


Always, <r 2 oc t 



SP, despite the inclusion of arbitrarily strong encounters, 
can be explained by the fact that the system evolves to its 
steady state much faster than the rate at which Ae > 1 
events contribute to relaxation. This can be shown in a 
quantitative way by defining an effective, time-dependent 
diffusion timescale for the SP. 

On long timescales, the SP is dominated by its central 
Gaussian, and therefore its width, a 2 — tT 2 log [L (q 2 bt)] 
(Eq. I34p , can be used to estimate an instantaneous diffu- 
sion timescale, tg s (t) — E^ ain /T 2 log L(q 2 bt). Therefore, 
at time t, the evolution of the system up to this time 
can be described by the FP equation with an effective 
(time-averaged) diffusion time 



1 



tf (t) 



1 



-/-111S 



■dt' 



(38) 



The analytic solution of the FP equation (see Section 
12 . 2(1 . t r w 10i£;(-E m i n ) then suggests to associate the time 
for relaxation to steady state by the SP with the solution 
of the ansatz t r ~ 10t c £(t r ), 



(E min ) log ( 



lOe 



3/2-7E 



/Ae n 



(39) 

which corresponds to an effective Ae max = 
^/I0e 3/2 ~ 7E ~ 8. This implies that by time t r , 
the SP can be represented by an FP propagator with 
a Coulomb logarithm that is formally increased to 
~ log(8/Ae m ; n ) (Figure [5]), thereby providing a lower 
limit on the time to relaxation due to the inclusion of 
strong encounters, t r > t F p log(l/Ae min )/ log(8/Ae min ). 

This lower limit is to be compared to the time for re- 
laxation that can be formally defined for the EP through 
its diffusion coeffici ent, w here for the p — 1/4 cusp, 
A w 56/Ae m ; n (see Eq. IA16|) . Assuming that the EP has 
already converged to the diffusive limit, this then implies 
that tf p /t s r p « log(8/Ae min )/log(56/A £min ) < 1. This 
is in contradiction with tf p being the lower limit on the 
time to relaxation, which means that the EP does not 
reach the diffusive limit by time t ri and therefore can- 
not be fully described by an FP equation. Nevertheless, 
the FP equation remains a reasonable approximation for 
describing relaxation around a MBH, since the discrep- 
ancy in timescales is not very large. For example, for 



(Section GO]), 0.76 
10 3 - 10 10 A'U 



< t PP /t b ; p < 0.92 in 



Ac m j n Q 

the range Q 

4. MEASURING ENERGY RELAXATION 

4.1. Overview of the best fit procedure 

The propagator W_B i j(A t £', t) (Eq. GUJ) is determined 
by the shape of the stellar cusp thr ough the known func- 
tions Y 1 {E, J) and T 2 {E, J) (Eqs. [A6j [A7J), and by the 




Fig. 5. — Time evolution of the perturbation's N p (x > 1) /N* 
for initial perturbation at xq = 10 5 as a function of time. The 
system is integrated by MC, with individual energy steps drawn 
from the FP propagator (triangles), the isotropic WP WE(AtE) 
(crosses), and the isotropic SP W^(AtE) (asterisks). The ana- 
lytically derived N p (x > 1) /TV* (Eq. I19H is shown for co mpa rison 
(solid line). The analytically derived N p (x > 1) /TV* (Eq. [T^J with 
e 3 / 2_ ' TE \/T0 is shown for comparison (dashed line). 



Ae ef f 

max 



parameters Ae m i n and Ae r 



However, on the short 



timescales accessible through our iV-body simulations, 
the propagator is fully determined by Ti, T 2 and the pa- 
rameter Ae m i n only (see Eq. [34| . In practice, Ti(E,J) 
has a vanishing contribution to energy relaxation on 
short timescales, so that the propagator can be approx- 
imated as a functional of T 2 with one free parameter, 

Ac m i n . 

As shown below, careful analysis of the A^-body results 
reveals that the propagator model provides an excellent 
description of the data (Figures [7] and [S]) . We thereby 
validate the functional form of T2, and determine the 
value of Ae m ; n . In principle, if it were also possible to 
validate Ti and determine Ae max , then the full shape 
and range of the transition probability (Eq. |2~4"|) could be 
determined, and the DCs T^ 2 = Fl 2 (E, J) log A could be 
calculated reliably. As it is, we have to rely on the fact 
that both Ti and T 2 were derived in the same theoretical 
framework, and therefore assume that the excellent fit for 
T 2 also implies a validation of T\. Furthermore, since in 
the small deflection angle limit Ae max ~ 0(1), setting 
Ae max — 1 is expected to have only a small effect on 
the numerical value of log A = log (Ae max /Ae m i n ). We 
find that with this empirically-motivated calibration, the 
numerical values of the DCs lie within < 20% of their 
theoretically predicted values. 
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4.2. Analysis of the N-body results 

To test and calibrate the analytic formulation of en- 
ergy relaxation, we carried out a suite of Newtonian 
A^-body simulations. Each simulation consisted of 10 4 
stars of equal mass M*, and a central massive object 
of mass M, = 10 6 M*. The initial semi- major axes of 
the stars were drawn from a p(a)da oc a 2 ^ 7 da distri- 
bution for 7 = 1.0, 1.25, 1.5 and 1.75, with random 
phases and isotropic orientations. The initial eccen- 
tricities were drawn from an isotropic p{e)de — 2ede 
distribution. In the Keplerian limit, these models ap- 
proximately correspond to a power-law density cusp, 
n(r) oc r~ 7 . The simulations were carried out with 
the parallel Nbody6++ code (implementing Hermite in- 
tegration with hierarchical block time steps, the Ahmad- 
Cohen neighbor scheme for efficient calculation of ac- 
celerations, and the Kustaanheimo-Stiefel regularization 
method for dealing with close few-body sub -systems) 
([Spurzem et al.|[200lHKhalisi fe Spurzem)l2006l) . 

Measuring energy relaxation in short-timescale N- 
body simulations is challenging, since the distribution 
of energy jumps is heavy-tailed, but the large energy 
exchange cutoff is not reached on these timescales (Sec- 
tion ^. 2p . Consequently, commonly used measures of evo- 
lution, such as the rms of the energy change, are strongly 
affected by rare strong exchanges and do not offer a ro- 
bust characterization of the evolving energy distribution. 
After some experimentation, we adopted a statistically 
robust method, described below, whose key idea is to 
estimate the width of the central Gaussian component 
of the energy distribution as the limit of fractional mo- 
ments. 

Snapshots of the A^-body simulation results are 
recorded at equally spaced times {ti = i<5i}i=o,i,2,...- 
For each star n of the total N, the orbital energy at 
ti, E™, is calculated and recorded. The stars are as- 
signed to logarithmically-spaced energy bins by their ini- 
tial orbital energy. For each star, and for each time- 
lag {A.tk = k5t}ts=i,2,3,..., the relative energy changes 
are collected on non-overlapping, consecutive intervals: 
{Ae£} = {(£«_ k - Ef )/E?} i= o, k ,2k,3k,..., thereby avoid- 
ing the complications of statistical inter-dependencies 
and correlations. The measured relative energy differ- 
ences, {Ae£}, are drawn from the heavy-tailed parent 
distribution W^e^j^Ae^, Atk), whose characteristic pa- 
rameters need to be estimated from the numerical data. 
We use the fact that the central part of W(Ae, t) is Gaus- 
sian to a good approximation with a scale parameter 
a\ E {t) = tT%{t) (Eq. [34]). Since most of the data lies 
within the central Gaussian (Figure HJ , its width can be 
robustly measured. As shown in Appendix [Bj the scale 
parameter can b e robustly estimated by fractional mo- 
ments 8 < 1 (e.g. lZumofen fc Kla fter 1994), which give a 
larger weight to the central part of the distribution and 
a smaller one to the diverging tails, 
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(40) 



where T(x) is the Gamma function and we assumed that 
WEiji (Ae, t) is symmetric since in our case, tT\ <C r 2 
(negligible drift). In particular, this relation also holds 
for a pure Gaussian (for any 8). We verified that the 
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Fig. 6. — The evolution of the energy scatter <r^ e bin( T ) (^1- 
1421 (blue circles) in iV-body simulations for several energy bins 
(displayed shifted on the y-axis for convenience) . On time-lags r > 
1 the scatter grows as r log r (red dashed lines) , as expected from 
our analysis. This is somewhat faster than the oc r rise expected 
for normal diffusion (black dashed lines). 

results converge rapidly as 8 — > 0; the results below are 
estimated with 8 = 10~ 3 . 

Over the short timescales of our simulations, E and 
J for each star can be approximated as fixed at their 
initial values E*, J*. Using Eq. [40l we define for each 
star a dimensionless scale parameter, which we measure 
directly from the data, 



^Ae,* = ~ = hill 



(|Ae a 



\S\2/S 



2n-V s T (Y 
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. (41) 

where Ai^ = tP*, P* = P(E+) and we use the minimal 
lag r = 1 in order to obtain maximal statistics (Figure 

We choose the energy bins to be wide enough so that 
the bin average over J is close enough to the isotropic 
average. Using Eq. (|4"Tj) . we define the bin averaged 
dimensionless scale parameter 



CT Ae,bi 



r%(E*,j*,TP*) T p* 

El 

r 2 (£ttn)TPbi> 



EL 



V 2Ae 2 . El 

\ mm bn 



where Ej,in and Pbin are the bin averaged values. The sec- 
ond approximate equality relates the measured quantity 
to the theoretical model (Eq. |3T>)) . under the assump- 
tion that the weak dependence of logL(a;) on x allows 
(xlogL(x)) bin to be approximated by (x) hin log L ((x) bin ) 
. Finally, we calculate T 2 from Eq. (|D2[) , and fit a\ t bin 
to the A^-body data to obtain the best fit estimate of 
Ae mm and a goodness of fit score for the model (Figure 
0). 

Our estimated best fit values for Ae min , averaged over 
~ 10 simulations each for various values of the power-law 
cusp, are presented in Table [51 We also list the relative 
errors resulting from the simple approximation A = Q, 
assuming Ae max = 1 . Figure [6] shows the growth of the 
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Fig. 7. — The energy diffusion rate as measured by 
°"L,bin ( T = U p ( E min)/P ( E bin) as a function of E bin /E min for 
different cusps. The analytical prediction (solid lines) is fitted to 
the measured values in the JV -bo dy simulations (dots), using Ae m ; n 
as a free parameter (see Eq. 1421 . 
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Fig. 8. — The distribution of energy jumps (gray crosses) mea- 
sured in the JV-body simulations, at one energy bin. R ed strip: 
exact numerical evaluation of the full propagator (Eq. 1291 , for the 
range of stellar angular momenta in the bin. Blue strip: the central 
Gaussian of the full propagato r (fi rst order term in Eq. 1341 1. Green 
strip: the FP propagator (Eq. 1321 . 



width of the central Gaussian, a\ e bin as function of the 
time-lag r. Pure diffusion would lead to a oc r behavior 
at t > 1 . However, the clear trend of a steeper than lin- 
ear growth ~ rlogr, a signature of anomalous diffusion. 

Finally, to verify our analysis we plot the histogram of 
the orbital energy jumps Ae^ in a given energy bin and 
compare it to the bin-averaged propagator 

W Ebin (Ae,t)~(W E ^(Ae,t)) Eun , (43) 

where we evaluated We+j* (Ae, t) nu mer ically for each 
star using the fitted AE^n, and Eq. (| A7|) . An example 
is shown in Figure[5] For comparison we also show the ex- 
pected bin-averaged FP propagator (W^ P j (Ae, t)) and 
the bin- averaged central Gaussian of the short timescale 
propagator (Eq. |3"4")) . 



5. APPLICATION EXAMPLE: STELLAR ORBITS 
AROUND THE GALACTIC MBH 

The dynamical effects of stellar encounters between 
stars in the G C may be detectable by next-generatio n 
telescopes fe.g. lWeinberg et al.ll2005tlBartko et ai1l2009D . 
and be used to probe the existence of the hypothesized 
stellar cusp there. Since the relaxation time in the in- 
ner GC is many orders of magnitude longer than the 
O(10 yr) observational timescale, such encounters will be 
in the anomalous diffusion regime. 

The question of the existence of the dark cusp near 
the MBH is closely tied to the puzzling paucity of the 
old luminous giants there, whose absence is in contradic- 
tion to the theoretical expectation for an old relaxed sys- 
tem (IBahcall fc Wolil 119761 \Wf% I Alexander fc Hopmanl 
120091) . A possible explanation is that a relaxed cusp 
does exist there, but for some reason, for example the 
selective destru ction of extended giant envelopes (e.g. 
IAlexander]|1999f ). the old giants do not trace the overall 
population. In that case the cusp is dark, that is, com- 
posed mainly of compact remnants and faint low-mass 
stars, and is expected to be dominated by ~ 10 M® stel- 
lar mass black holes that have accumulated by the pro- 
cess of strong mass segregat ion over cosmological times 
(j Alexander fc Hopmanl 120091 ) . 

Changes in orbital energy can be expressed as changes 
in the radial orbital period (equivalently, a phase drift). 
Assuming that the dominant cause of such changes is 2- 
body interactions, or that it is possible to model and 
subtract all other competing effects, then one way of 
detecting or constraining the dark cusp is through the 
stochastic orbital phase drift in the few observed lumi- 
nous stars very close to the MBH. 

Figure[5]shows the probability of a period change dP/P 
after 10 years for a test star with an orbital period of one 
year (semi-major axis of 0.7 mpc), as a result of inter- 
acting with a relaxed main sequence cusp (p(r) ~ r -7 / 4 ) 
and with a mass-segregated dark cusp (p(r) ~ r 3 / 2 +P) 
of stellar remnants. Figure [TU] shows that a change of 
dP/P = 0.0015 over 10 years due to a mass segregated 
dark cusp, measurable by next-generation telescopes, has 
a non- negligible probability of sa 3.3% per star. It should 
be noted that assuming only the effects of strong single 
encounters leads to an under- estimate of the probabil- 
ity (1.5%) of such an event, while the FP probability 
of 10.5% is an over-estimate, because the event lies at 
the relatively under-populated intermediate wings of the 
propagator (Figure The naive FP estimate could 

therefore lead to erroneous conclusions about the ab- 
sence of a cusp in the case of a null result. For exam- 
ple, if an event of dP/P = 0.0015 is detected, the FP 
{l,2,3,4,5}er confidence levels on the minimal number 
of stars in the cusp, min ncr Npp, will differ from the true 
ones (due to anomalous diffusion), min n(T N, by factors 
min„ CT iV> F /min n<r iV = {0.52, 0.55, 1.7, 38, 2.7 x 10 3 }, 
respectively. Thus, the erroneous assumption of normal 
diffusion may result in either over- or under-estimates of 
the dark cusp density, depending on the chosen confi- 
dence level. 

6. DISCUSSION AND SUMMARY 

The stochastic evolution of a stellar system is described 
at the most basic level by the transition probability, 
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TABLE 2 
Measured Ae 



7 




A/Q C 1 ' 2 ) 


log A/ logQ^ 1 - 2 ) 


X^ r H ( A(E min) 




1.0 


5.3 ±0.5 


0.19 ±0.02 


0.88 ±0.01 


1.2 ±0.3 


8 


1.25 


6.7 ± 0.7 


0.15 ±0.02 


0.86 ±0.01 


1.3 ±0.2 


12 


1.5 


9.2 ± 1.3 


0.11 ± 0.02 


0.84 ±0.01 


1.5 ±0.5 


9 


1.75 


13.7 ± 3.6 


0.07 ± 0.02 


0.81 ± 0.02 


2.0 ±0.6 


9 



1 For Q = 10 6 , best fit values for JV-body results f Section IP! , 

2 A = Ae max /Ae min , where Ae max = 1 




dP/P 

Fig. 9. — Probability for a change in the orbital period of a 
star with a 1 year period, over an observation time of 10 years, 
for several background population models. The probability for an 
accumulated period change due to multi ple energy changes calcu- 
lated by the short time propagator (Eg. 129( 1 (solid lines), is larger 
than the probability for this change due to a single energy scatter 
(dash-dotted lines). For comparison, we show the accumulated pe- 
riod change probability calculated by the FP propagator (dashed 
line). The background population is modeled by a cusp of IMq 
stars with a slope of 7 = 1.75 (green) with a total mass of IO^Mq 
within 1 pc, and by a cusp of IOMq stellar remnants with a slope 
of 7 = 1.9 (blue) and a total mass of 10 5 Mq within 1 pc. 
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Fig. 10. — Visualization of the orbital perturbation of a one year 
Keplerian orbit with an eccentricity of 0.9 to a level of AP/P = 
0.0015 (blue), 0.003 (green) and 0.006 (red), sampled on equal 
times, with respect to the original period (black). The maximal 
spatial shifts are 0.415 mas, 0.824 mas and 1.61 mas. 



which expresses all the physical processes considered, and 
all the approximations assumed. The accumulated effect 
of multiple scattering events over a time that is short 
enough for the system to remain almost unchanged, is 
constructed by integration over the transition probabil- 
ity, and is described by the propagator. This general 
description of the short- or medium-term evolution can 
lead to a wide range of stochastic processes. Only in 
the specific case when the propagator asymptotes to a 
Gaussian, which is fully defined by its two first moments 
(the first and second DCs), does the evolution progress 
as cx y/t and the process becomes normal diffusion, which 
is described by the FP equation. However, it is not at 
all guaranteed that this convenient and useful limit pro- 
vides a valid description of the accumulated effects of 
the transition probability — this has to be proved case by 
case. There are in fact many examples of anomalous 
diffusive processes, whose origins can be traced to diver- 
gences of the moments of the transitio n probability, such 
as th e super-diffusive Levy flights (e.g. lMetzler fc Klafterl 
l2000h 

The FP equation has been applied extensively in past 
studies of nuclear dynamics near a MBH. However, the 
simplifications, assumptions and free parameters (many 
necessitated by the diverging, un-screened nature of the 
gravitational force), which are needed to justify and use 
this method, remained largely untested. This called for a 
more rigorous validation of the basic underlying assump- 
tions. This need was further motivated by some con- 
flicting theoretical and numerical results on the emerg- 
ing time scales of diffusion and relaxation in galactic nu- 
clei (e .g. iBahcall fc Woilll976t iPreto et alJl2004t iMerrittJ 
l2010t Madiga n et al.l 120111) . These impact multiple is- 
sues of interest: the rates of tidal disruptions, the rates 
of GW emission events from inspiralling compact objects, 
and the stellar distributions very close to MBHs, and in 
particular those of compact remnants and stars on rel- 
ativistic orbits. Furthermore, the one system where the 
stellar distribution around an MBH can be directly ob- 
served, the GC, does not follow the FP predictions for an 
isolated old population, in contrast to prior expectations 
dBuchholz et al.l l2009l IDo et al.ll2009t iBartko et alllMToL 
lMerrittj|2010f ). 

Past studies of relaxation theory usu ally bypassed th e 
use of the transition probability (but see lGoodmadll983f) . 
They assumed from the outset that the diffusion coef- 
ficients exist and only the first two are important, and 
calculated them by adopting cutoffs (the small angle and 
hyperbolic orbits approximations) for the exchanged en- 
ergy in 2-body interactions. However, this procedure 
limits the possibility of describing the behavior on short- 
timescales, of verifying that relaxation does in fact con- 
verge to normal diffusion, and of estimating the contri- 
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bution of strong encounters. 

We focused here on the short-timescale behavior of en- 
ergy relaxation, which can be directly probed by high- 
accuracy iV-body simulations, and is easier to describe 
analytically in terms of small deviations from an almost 
stationary DF. We derived a rigorous analytical descrip- 
tion of relaxation on short-timescales, which we validated 
by fitting to numerical results using a robust statistical 
analysis method. This then made it possible to tie the 
short-timescale evolution, calibrated by the iV-body ex- 
periments, to the asymptotic FP one, thereby fixing the 
absolute timescales for diffusion and relaxation around a 
MBH to within 0(10%). It should be emphasized that 
the convergence of the propagator to the central limit, 
and the return of the stellar system to steady state, are 
in principle unrelated processes, which proceed at unre- 
lated rates. There is therefore no a-priori guarantee that 
the FP treatment is a good approximation for describ- 
ing relaxation. However, we showed that in the weak 
limit, the evolution of the propagator converges to the FP 
description long before the system relaxes to its steady 
state. We also showed that the inclusion of strong col- 
lisions shortens the convergence time by less than a fac- 
tor of 4/3, because the system reaches steady state be- 
fore these rare collisions have time to play a substantial 
role. Together these results validate the use of the FP 
approximation on the intermediate and long timescales, 
and prove that a perturbed BW cusp will regain is steady 
state after t r ~ lOi^, where fg is the diffusion time at the 



outer boundary. We conclude that stellar cusps around 
lower-mass MBHs (M. < 10 7 Af Q ), which have evolved 
passively over a Hubble time, should be dynamically re- 
laxed. Anomalous diffusion can affect processes observed 
on short timescales. As an example, we showed that if 
future observation of stars orbiting the Galactic MBH 
detect orbital energy perturbations due to the dark cusp 
of stellar remnants that is expected to be there, estimates 
of the cusp density will depend critically on accounting 
for the anomalous nature of the diffusion. 

Several open question remain. In this study we 
restricted ourselves to energy relaxation in a single mass 
population. A mass spectrum will generally change 
the steady state and s horten t he relaxation time (e.g. 
Alexander fc Hopmanl l2009f . lAmaro-Seoane fc Pretol 
20101 ). It is likely that the convergence of the propaga- 
tors will depend on mass, and the use of the multi-mass 
FP equation will then have to be re-examined and 
justified. The evolution of angular momentum, which 
is crucial for replenishing stars destroyed by the MBH, 
and thus competes with energy relaxation, was not 
addressed here at all. The same considerations of 
anomalous diffusion should in principle apply also to 
angular momentum. This requires further study. 
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APPENDIX 



A. THE ENERGY TRANSITION PROBABILITY 



We derive here explicit expressions for the energy tran sition probability for several specific cases. The energy 
transition probability is generally given bv lGoodmanl (119831 Eg. AllR 



K r (E,E') = 



2^G 2 m 2 J7j\„ dE f f (E f ) (§ W 2 - AAE) v + r|«+ A£ dE f f (E f ) (§w* - AAe) v' f AE < 



v\AE\ 



3 

f^dE f f(E f )(§v» 



- 4AE) v> + J^ AE dE f f (E f ) + 4AE) v f AE > 



(Al) 



where E' = E + AE, v f = y/2(<j)-Ef), v' f = v / 2(cj) - Ef + AE), v = \/2(<j) — E) and v' = y/2(<j)-E'). Note that 
here E is the positively defined orbital energy and <j> (r) is positively defined potential. Thus Eq. (|A1[) reads 



K r , E {AE). 



32 ir 2 G 2 m 2 , 
IT 



3 \AE t 

' A 3 AE 



4 ip-E 
ip{r)+AE 
E+AE 
1 AE 



)j^ E dE f f(E f ) 
d Eff(Ef)^) 



3/2 



1 AE 



4 tp-Ef 



1 - 



AE 



rE+AE 



4 ip-E J Y V-E JE 

£l E dEff(E f )(^ 



dEjf(Ej) 



3/2 



3 AE 



4 ip-E s 



AE 



AE < 



AE > 



(A2) 



2 We correct here a couple of printing errors in Eq. All of 
I Goodman] 111983 ). The correct prefactor should have a tt 2 term, 



and not n (cf Eq. A7 there). The first term of the second line 
should have v' terms and not v (cf Eq. A12 there). 
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Assuming \AE\ < E and expanding Eq. (IA2[) in AE, we obtain 
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Taking the orbit average Kj_e{AE) = 2P 1 f™ K r ^E{AE)dr /v r , we obtain 



K E j(AE) = 

' Jy ' 2 \AE\ 
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Assuming that the DF is isotropic in velocity space, by averaging over angular momentum 

{X) lso = ^ / C Xdj2 = 4 / Xr 2 v r dv r , 

J c Jo J c JO 

we obtain 



T 1 (E) = 
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where g(s) = f{sE)/f{E). 

An isotropic averaged transition probability can be directly obtained from Eq. (|A2j) without assuming 
limit, 

K E {AE) = ^f dv r r —r 2 v r K r , E (AE) = -^-~H (Ac) , 



EJc Jo 



(A3) 

(A4) 
(A5) 
(A6) 

(A7) 

(A8) 

(A9) 

(A10) 

the weak 
(All) 



where K E = 3 -ir L G 2 m 2 f (E) E 2 and 



(1 - | Ac) g (s + Ae) ds + (a - f Ac) ^g (a + Ac) Ac < 
SL 9 (s + Ac) ds + /r (s + 1 Ac) jj^g (a + Ac) Ac > 
Since isotropic averaging includes J — > where Ae diverges, the DCs are given by 



If (E) = K E (E) Qf° [H (Ac) - H (-Ac)] 
If (E) = K (E) QH [if (Ac) + (-Ac)] 



dAe 
Ac 2 " 

dAe 
Ae 



After integrating by parts, taking the limit Ae m i n — > and using the relations I^ (i?) = K-eH (0) and 
K-eH' (0), the DCs can be written as 



(A12) 

(A13) 
(A14) 
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Note that for a power law distribution, g (s) — s p , Eq. (|A12j) reads 

f ^ (- A £ f- 1/2 (5 - 2p) r (| - p) r (1 + p) Ae < -1 

H (Ae) = | 1 ja+A.) 1 ^ + (5 _ 2 p)(-Ae)-*+P / 3 (-Ae, | -p, 1 +p)] -1< Ae < (A17) 

i^ 4(3 " 2p) ^ ( )r 4p)1A£ (^ + Ae) p ~ 3 / 2 A, > 
where /? (x; a, 6) is the incomplete beta function. 

B. EVALUATION OF THE STRONG PROPAGATOR 

We derive here a perturbative expression for the SP in the limit AE max — > oo. In this limit, which is the relevant 
one for evolution on short timescales, the problem is closely related to the problem of Pareto sums with index a = 2. 

Consider a star with binding energy E and angular momentum J. The orbit-averaged rate with which the test star 
changes its energy is given by the transition probability (Eq. I24[) 



When the timescale of interest is much shorter than the energy relaxation time, Fi and T 2 can be approximated as 
constant in E and J and the maximal energy change in a single scattering event is much smaller than the large energy 
exchange cutoff AE -C AE^*. Therefore A£" max can be taken to the limit AE max — > oo. In that limit the Fourier 
transform of Ke.j(AE) is given by 
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JAE mir , 
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= q 2b {1 + m (2 - 2 lE - log(fc 2 A^ in )) (kAE min ) 

- \ (3 - 2 lE - log(fc 2 Ai£ . J) (kAE min ) 2 + O (fc 3 A^ in ) | , (B2) 

where ?7Ae min = ET^E, J)/T 2 (E, J)Ae min < 1, Ae min = AE miu /E, q 2b = J K EJ (AE)dAE and j E is the Euler's 
gamma constant. The energy propagator can therefore be directly obtained from Eq. (|29p 



~ tkAE dk 



W(AE,t) = ±-l 

exp {77Ae min i (2 - 2j E - \og(k 2 AE 2 nin )) kAE min q 2b t 

- \ (3 - 2 lE - log(fc 2 AE 2 min )) (kAE miD ) 2 q 2b t^ . (B3) 



Define r = <72&£, c(r) = ^/ 2r log L(r) Ai?^ in and /x(r) = 2?7Ae m i n r logL(r)A£ , min , where L(r) is defined by 

2 log L (r) - log(2 log L (r)) = 3 - 2 lE + log(r) , (B4) 

and 

L (t) = exp f -^-) \ M V^exp (l + l0 f ^ 1 , (B5) 



log (e 3 - 2 ~t E T) 

where W—i(t) is the Lambert W function, and L is an increasing function of time. By a change of variables to 
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x = (AE - (x(t))/(t(t) and uj = fccr(r), Eq. (|B3| reads 

,2l„„,.,2 



V21oiW (l + 10gCj2) ) W 



i r , ,.-r. , - 2 io§- 2 if?AeD 



2W-oo I 41ogL(r) 

1 = e ~^ + 7~i — ~FTT fe^) + 2^Ae min ^4^ c 2 (x) ) , (B6) 



where 



and 



V^F 41ogL(r) V CW *(t) 

W c x (a;)= T- 1 {e-^a; 2 log(w 2 )} = 7/= e ~^ { 2 + - X ) [^(a; 2 ) + j E + log 2] } - serf f J|0 (B7) 

V^ 2 (a:)=^- 1 {- l e-^ W (l + logV))} = -j=e-^x( lE ~l + \og2 + F 2 (x 2 )) - erf f-?J=) , (B8) 



where -F 2 (a:) = a^i ^CE -*-) 3/2, 2,x/2) and 2-F2 is a generalized Hypergeometric function. 
Integrating Eqs. (|B7[) . (|B8j) to obtain 



and 



FWl{x) = J Wl{x')dx' =F~ 1 \iue-^\og{u 2 )} = erf - -^=e~^ a; ( 7£; + log2 + ^(a; 2 )) , (B9) 

FW/ 2 (x) = / W?(x')dx / =^ 1 {e- i ^(l + log(w a ))} = — i=e~^ [ 7B - 1 + log2 + F 2 (x 2 )] (BIO) 

J-oo ' V27T 



and therefore 

1 

1 + erf 

' 2 

By Defining 



( FWl(x) + 2r 1 Ae min ^lFW 2 (x) 



41o gj L(r) V ct(t)' 



rf ia (B,j,t) = r 1)a (B,j)iogL(ga6(s,j)t) , (Bii) 

we can write (i (t) = tTf(E, J, t) and a 2 (t) = tT%(E, J, t) and Eq. ([B6|) reads 

y/2%tT%{E,J,t) 4 log L(g 26 t) y/tT%(E,J,t) c \ y/tT%{E,J,t) J 



1 1 /r^^/Ag-^^\ (m2) 



V41ogi(g2fct) y/tT%(E,J,t) V £ 2 c ^ y/tT%{E,J,t) 

This is a modification to the impulse response of the FP propagator (Eq. (|32p) that is a Gaussian with a variance of 
£T2 log A. If we assume that Ai? max ~ E and AE m i n ~ A _1 _E, then the time £a for which the Gaussian in Eq. (|B12|) 
has the same width of the FP Gaussian is 

t A = e 2 ^~ 3 — f— « t B /6. (B13) 

r 2 log a 

To the first order, the asymptotic expansions of and W 2 are cc/|a;| 3 , 2/|a;| 3 respectively, and thus all the moments 
diverge. The mean and variance are therefore not useful estimators. However, the fractional moments (i.e. [i$ = 
(\AE\ 6 ) where S < 1) do converge, which suggest that these moments can be used as estimators. Using Eq. (IB12[) the 
fractional moments is defined as 



fie(T) = J dAE\AE\ 5 W B ,j(AE,t) 
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In the limit 5 < 1 and/i 2 /cr 2 < 1 



2/5 



-2a 2 





"i + <r 








) 




2 





2/<5 



(B15) 



and therefore cr 2 can be estimated using 



cr 2 = lim 



1 



(AE) 



s-o 2 \ r (Y; 



2/5 



(B16) 



C. TIME-DEPENDENT SOLUTIONS OF THE FP EQUATION FOR A POWER-LAW CUSP 

We solve Here analytically the FP equation describing the evolution of a small perturbation on a fixed background 
distribution of a power-law cusp around a MBH. Consider a small perturbation N p (E) on a fixed background distri- 
bution Nb(E) oc E p ~ 5 / 2 . Assuming an infinite cusp, the DCs are given by T^(E) oc EP+ 1 and T$(E) oc E'P+ 2 . Define 
x = E/Eh, where Eh is some reference energy, for example Eh = a 2 where oq is the velocity dispersion far from the 
MBH. 

The FP equation for the small perturbation N p is 



9N P (x, t) _ l_d_ p+2 
dr 



(Cl) 



where 77 = E h T^{E h )/T^{E h ) = -(1 - 4p)(l - 2p)/(4(3 - 2p)) and r - t/t| where t% = E 2 /(T 2 (E h ) log A). The 
steady state solution (assuming p + 2 — 2i] > 0) is N p yD (E) oc E 2r, ~ p ~ 2 . Assuming particles cannot leave the system 
and are restricted to energies E > Eh, the steady state is given by N£° (x) = (p + 1 — 277) cc 2?)_p_2 . 

The time dependent solution can be obtained by the eigenfunction method (e.g. lGardinerl l2004f). Let N p (x,t) = 
N p >D (x)q(x,T). We are looking solutions of the form N p (x,t) — P\(x)e~ Xl ~ and q(x, r) = Q\{x)e~ Xr , which satisfy 



1 



,P+2 



and 



F p (x) 



o2 a 

— Q A (a;) + vx p +1 —Q x (x) = -XQ 



p+l-2rj 



dx 



x 2r >e~ XT - 



.(*), 



dx 



Q(x). 



For p > the solutions of Eqs. f([C2]l, (fC3)0 is of the form 



(C2) 
(C3) 

(C4) 



P A (a;) = C Aa ;^- 3 / 2 J Q 



2V2A 



pa; 



p/2 



Qx (x) - 



p + 1 - 27/ 



2V2A 



px 



p/2 



(C5) 
(C6) 



and 



F„ 



{x) = Cx2V2\ x _ {1+p+2h) /2 Ja+i 



2\/2A 



px 



p/2 



(C7) 



where a = (1 — 2r\)jp, and J«(x) is the Bessel function of the first kind, and where we assumed F p (x) = as x — > 00. 
The assumption that stars cannot leave the system and are restricted to energies E > Eh is expressed by a reflecting 
boundary at Eh, i.e. F p (x = 1) = 0. The eigenvalues are therefore 



A,: 



o Ja-\-l,n ) 



(C8) 



where j'„ in is the n-th zero of J v {x). 
By choosing Ca = -p, the eigenfunctions Q rl = Qa„ and P„ = Pa„ satisfy the bi-orthogonal relations 



PnOzOQmOr) 



(jo:+l,n) (ja+l,m) Jo 



sJa. (ja-j-l,n s)Ja(j 



a+l,m 



s)ds 



(C9) 
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Any solution can therefore be written in the form 

N p (x,t) =^A n P n ( 



x)e 



-A„t 



(CIO) 



where A n = J^ 1 Q n N p (x,0). 

The time dependent distributions resulting for initial conditions of N p (x, 0) = S(x — xq) ) for example, are given by 



N p (x,t\x ,0) = N p x (x) 



1 + Qn{x)Qn(xo)e-^ j *+i." T 



where 



Qn (x) 



Ja(ja+l,n) V P + 1 ~ 2? ? 



-J, 



(j Q+ i,„*- p/2 ) 



(Cll) 



(C12) 



A different situation is obtained when the test stars can leave the system. This can be approximated by imposing a 
boundary condition N p (x — 1, t) — 0. In that case A = p 2 j& „/8 and 



Q(x) = 



p+l - 2r) J a + 1 (j a .n) 

The response for <5(a; — xq) initial conditions is then given by 



iV p (x,t| a; o,0)=iV p 00 (x)^Q„(a ; )Q„(xo: 



e * 



(C13) 



(C14) 



where the number of test stars with relative energy larger than x at time r is given by 



N p (>x,r) =2(x/x ) 1/2 - v x-^ 2 Y, 



I _li f? x^/ 2 ! 7 (i x~ p/2 \ 



p -2 



(C15) 



D. ENERGY DCS FOR A FINITE ISOTROPIC POWER-LAW CUSP 

Cusps around MBHs, whether real or in A-body simulations, contain a finite number of stars, and are therefore 
finite in extent. As we show here, these edges can substantially increase the relaxation timescale. 

Consider a system with N stars of mass m each, with a semi-major axis distribution p(a)da oc a 2 ~ 1 da and eccentricity 
distribution p(e)de = 2ede. The number of stars with energy larger than E is N(> E) = N tot (E min / E) 3 ~^ where N tot 
is the total number of stars up to energy E m { n and 7 = 3/2 + p. Since the number of stars in the system is finite, the 

DF is t runc at ed at some maximal energy E max w ^tat 
Eqs. {59), (|AT0| imply (for 0.5 < 7 < 2) 



NW J, E min , where N(< £ max ) = 1. 



If (£) = log A-^—^ N(> E)S E1 {E, E min , E max ) 



P{E) Q 2 



and 



where 



r£(£0=iogA 



E 2 c 2 
P{E) Q 2 



Cl 



cS = 16 



N{> E)SE22{E,E min ,E n 



4(1 - Ap) 



l+p ' 
3-2p 



(l+p)(l-2p) ' 



(Dl) 



(D2) 



(D3) 



(D4) 



and Sei, Se2 a- re the correction functions due to the edges, whose parameters are p and N tot (via E max / E nlul = 



N, 



l/(3/2-p) 



2+2p 
l-4p 



3/2-p 



Ssi(£) = < 



1 - 



3-2p ( E mi 
l-4p V E 



3-2p J 1 
l-4p ' 1 



P+ 1 , 2+2p / g \ 

"T" l-4p ^Bmax^ 

(fc) 1 " 



3/2-p 



3/2-p 



-E < -^min 

-^min ^> E ^> -£/ max 



(D5) 
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Fig. 11. — Energy relaxation as function of energy, for a cusp (N(> E) o c E 3 ~~ l ) with Ntot = 10 4 stars and Q = 10 6 , fo r different power- 
law sl opes. The solid lines are the relaxation rates calculated by Eq. ID2t . the crosses are the rates calculated using the Colin & Kulsrud 
( 1978) integrals, the dashed lines are the rates for infinite cusps. 



and 



Se22{E) — < 



2+2p 
3 

1 - 



(A) 



1/2-p 



l-2p ( E n 



l-2p (E mssL \P+ 1 



1 - 



(fc) 



_ 2+2p 



1 



grain 
£max 



1/2-P 



p+1 



1/2-p 



-^min ^> E ^> E max 
E > E max 



(D6) 



3 V E 

Figure QT] shows that the truncation of the DF in a finite cusp leads to slower relaxation 
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